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In the presence of a magnetic field and an external periodic potential, the Landau level spectrum of 
a two-dimensional electron gas exhibits a fractal pattern in the energy spectrum which is described 
as the Hofstadter’s butterfly. In this work, we develop a Hartree-Fock theory to deal with the 
electron-electron interaction in the Hofstadter’s butterfly state in a finite-size graphene with periodic 
boundary conditions, in which we include both spin and valley degrees of freedom. We then treat the 
butterfly state as an electron crystal so that we could obtain the order parameters of the crystal in the 
momentum space and also in an infinite sample. The excitation gaps obtained in the infinite sample 
is comparable to those in the finite-size study, and agree with a recent experimental observation. 


In a strong perpendicular magnetic field, the energy 
spectrum of a two-dimensional electron gas (2DEG) splits 
into a series of Landau levels. If a periodic potential 
is also added to this system, the noninteracting energy 
spectrum then displays an intricate fractal pattern known 
as the Hoftstadter’s butterfly flj. In the quantum Hall 
effect regime, the Coulomb interaction between electrons 
plays an important role both in the ground state and in 
the low-energy excitations. Theoretically, the Coulomb 
interaction effects on the butterfly states were investi¬ 
gated in the Hartree or mean field approximation 
or with the electron correlations @, 16 ]. In the Hartree 
approximation, the electrons are classical particles with 
negative charge and repel each other. This approxima¬ 
tion is unable to deal with a spin/valley system such 
as graphene if the spin is not polarized, since the ex¬ 
change interaction is not taken into consideration. In 
the exact diagonalization scheme, the electron correla¬ 
tions are completely included, and hence more adaptable 
to the fractional quantum Hall effect Q, but then only 
the finite-size systems can be handled in this scheme. 
Here we consider the Hartree-Fock approximation (HFA) 
to deal with the Coulomb interaction and the spin/valley 
system in graphene. Further, we work in the momentum 
space to study an infinte system. 

After a moire pattern is fabricated by two misaligned 
honeycomb lattices by graphene and the boron nitride 
(BN) substate mm , the fractal band structure and the 
transport properties of the butterfly states were finally re¬ 
vealed in recent experiments |12j, |Tj] . Theoretical studies 
of the butterfly states in monolayer graphene and bilaver 
graphene have exhibited very rich physics asm. 

In contrast to the conventional 2DEG in a semicon¬ 
ductor, electrons in graphene must be described by the 
Dirac equation with an extra valley degree of freedom 
fl9U2il . In this case, the Hartree approximation may 
not be enough to describe the behavior of the electrons 
if we consider both the spin and the valley. Indeed, the 
excitation gap measured in a recent experiment fl3j can 
not be explained either in the noninteracting picture or in 
the Hartree approximation. Here, we develop a Hartree- 
Fock approximation to calculate the energy gaps for in¬ 
teger filling factors. We also employ a method involving 
the crystalline state to explain the experimental results. 


The Hofstadter states are not affected much by the 
geometry of the external potential. For simplicity and 
without loss of generality we introduce in the graphene 
Hamiltonian a square periodic potential [23j , V/ X t (£, y) = 
Vq [cos ( qox ) + cos (qoy)] , where Vo is the amplitude of 
the potential and qo = 2tt /a® with a period of the poten¬ 
tial a o- The Hamiltonian of the 2DEG in such a system 
is then given = H 0 + Vc + T4xt ■ The noninteracting 
Hamiltonian without the external potential is 


H 0 = v F 


0 

P x + ir]P y 


P x - ir\P y \ 

0 J 


(1) 


where r) = ±1 for K and K' valley respectively, v F is 
the Fermi velocity, and P = p + eA is the canonical 
momentum. We choose the Landau gauge to define the 
vector potential A = (0 ,Bx), and Vc is the electron- 
electron interaction. 

The period of the external potential is large enough 
experimentally mm (larger than lOnm), so that the 
valley mixing can be neglected. The energy bandwidth 
is narrow when Vq is not large. The Landau level (LL) 
mixing is also very weak and is not considered here. We 
first diagonalize the noninteracting Hamiltonian with the 
external potential, Hq = Hq + I4xt by using the basis 
{4>x}, where (f>x is the wave function of the n-th LL 
with valley-spin index er and the guiding center X [U- 
[22l ]. Then we could obtain AN^ eigenvectors: 



■i,x N ^x N . 



where i = 1,.. .AN,/, and 4 is the degeneracy of the 
Landau level in the finite sample. The coefficients cf x 


satisfy the normalization condition Sj=i H x, 

1. Note that i is also the index of the corresponding 
eigenenergies with ascending order and j is the guiding 
center index for the N$ states in a single LL. 


The Hartree-Fock Coulomb interaction shall be deter¬ 
mined self-consistently by the coefficients c? Xj - Let us 
consider the Hartree term V F and the Fock term V F sep- 
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arately. The Hartree term i2f, HH is 


(a',X'\V H \a,X) 


r 2 ' ___ N <j> 

Oa,a' c \ ' \ ' \ ' \ ^ 

V i G £ k,l 


1 

G£ 


*4!kii M k,i (— G) M x > >x (G), 


( 3 ) 


where n is the dielectric constant, the summation with a 
prime sums over all the states below the Fermi level, and 
the bar over the summation means the term with G = 0 
is excluded. The function M is defined by 


(G) 


&k,l+G y P i Gx ( k+ i) G 2 e 2 /4 
2lSn,o 6 


X 



J \n\ 



( 4 ) 


with a Laguerre polynomial L (we define L n< q = 0) and 
the magnetic length £ = ^JKfJeB). The Fock term is 


calculate the energy spectrum in an infinite sample, we 
could work in the momentum space. To be consistent 
with the finite-size study, we consider a crystal phase of 
the electron gas with the same geometry as the periodic 
potential. The lattice constant of this electron crystal 
a is given by a/1 = \j2'Kn c /v^ where n c is the electron 
number per crystal site. If the lattice constant of the 
electron crystal is identical to the period of the potential 
a = ao, then the electron number per site is given by 
n c = v la. 

The Hamiltonian of the 2DEG in the Hartree-Fock ap¬ 
proximation is written 

H = E a p a j(7 (q = 0) + 4^ext (q) Pa,a (q) (6) 

a C7,q 

+E E Uh (~ q )) (q) 

q a,a' 

~ EE Ux ( q ) (~ q )> (q) ■ 

q a,a' 


Nt 


(a',X'\V F \a,X) = 


EEE 


ntm ^ ^ ^ Gi 

V i G k.l 


xc i*k c i,i Mk ’ x (~ G ) M X',1 (G). 


( 5 ) 


For a finite sample the momentum vectors are discrete, 
G = (2n/ao) ( n x ,n y ), where n x ,n y are integers. 

Once we diagonalize the noninteracting Hamiltonian 
H 0 , we obtain a series of coefficients cf A -.. For the 
i- th eigenenergy Ei, the eigenvector is given by c° x .. 

The summation Y4 in the Hartree-Fock approximation 
contains i = 1 ...N s . The Filing factor is defined by 
v = Ng/N^ . We use these coeffeicients c? x . to com¬ 
pute the Hartree-Fock Coulomb interaction. Then the 
full Hamiltonian H is diagonalized and a new group of 
the coefficients d* Xj is obtained. By repeating this pro¬ 
cess the coefficients cf x and the energy spectrum can be 
evaluated self-consistently. The energy gap is obtained 
from A = E Xa+ i — E Xa . It is the gap between the high¬ 
est occupied state and the lowest unoccupied state. This 
gap may be observed in transport or capacitance mea¬ 
surements. 

The parameter a defines the units of the magnetic flux 

per unit cell of the periodic potential, a = $/$o, where 
<f) 0 is the magnetic flux quantum. Hence, a also describes 
the magnetic held if the periodic potential is fixed. More¬ 
over, the size of the sample is related to the period and 
a : = L x L y / (aag) , where L x and L y are the length 

of the sample in the x and y directions. In what fol¬ 
lows, we would like to study the energy gaps in different 
magnetic fields or a for a fixed Riling factor v. However, 
the size of the sample is not fixed ( L x and L y are not 
constant) when we study the system for a continuous a. 

Crystal phase in an infinite sample: If we calculate the 
energy spectrum for continuous a in a fixed-size sample, 
the energy gaps may not be reliable, since both the en¬ 
ergy spectrum and the gaps may be size-dependent. To 


The density matrix operator is 


?,a' (q) = E E e lq ^ X+X ^ 5 x,x-+ qy pcl x c a ^ xl , 

$ X,X ' 


( 7 ) 

where operators c a} x , 4 x are the annihilation and cre¬ 
ation operators of electrons in valley-spin a and the guid¬ 
ing center X. The Hartree and Fock interaction func¬ 
tions, Uh and U x , are defined by 


Uh( q) 
Ux (q) 


e 2 1 ■ 

[ M q v e 2 /2,- qy P/2 (<?)]' 


k£ , 


dp Mpy_ _Py_ (p/l) Jo (pq£) 


( 8 ) 

( 9 ) 


where Jq is the Bessel function. The external potential 
in such an electron crystal is 

Kxt (q) = yM, s< 2 /2 ,- qy p/2 (?) ■ (10) 

We define the Green’s function G a ,a’ {X, X', r) = 
~(rc a ,x (r)4, iX , (0)^ , where T is the time order oper¬ 
ator. Then at zero temperature, 

Ga,a' (q,T = 0") = {pa',a (q)) , (H) 

where G a ,a' (q,r) is the Fourier transform of 
G a ,a' (X,X',t). In the Matsubara frequency w n , 
the equation of motion of the Green’s function 


( ihuj n Ed) Gd,e (q,W n ) — hSd,eh q,0 (12) 

+ ]T Hext(q')e Jq ' Xq ' 2/2 G d ,e (q + q',w n ) 
q' 

+ E Uh ( q/ ) (/W (~ q/ )) e * q Xqr/2 Gd,e (q + q',Wn) 

a,q' 

- E U x(<l') (Pa,d (-q')} e iq ' xq ^ 2/2 G 0 ., e (q + q',w„), 

a,q' 
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where d, e are valley-spin indices, can be solved self- 
consistently (see [26] for details). Then all the elements of 
the density matrix can be obtained according to Eq. (ED- 

The density matrix completely describes the system 
with the Hamiltonian in Eq. ©• The energy spectrum 
is thus obtained by solving the equation of motion in 
Eq. (fT21) . When we estimate the energy gap, we need 
to subtract the energy of the lowest unoccupied state by 
the energy of the highest occupied state in the density of 
states (DOS). The relation between the DOS g and the 
retarded Green’s function is [27[ 

sM = [ G £<7 -> W + J0+)] . (13) 

<J 

The crystallized electron gas was studied in monolayer 
graphene [25| and graphene bilayer [27j without the ex¬ 
ternal potential. Wigner crystals and skyrmion crystals 
can be found in those systems at non-integer filling fac¬ 
tors. Moreover, the skyrmion crystals were found in bi¬ 
layer [28} and trilayer graphene [0 even for integer filling 
factors, due to the Dzyaloshinskii-Moriya (DM) interac¬ 
tion. We employ the same method to study the electron 
gas in the presence of the external potential. For a strong 
potential the 2DEG may be crystallized with the same 
geometry as the potential for integer filling factors even 
without the DM interaction. 

In a recent experiment 0 , there is about 1° misalign¬ 
ment between graphene and the BN substrate with the 
dielectric constant k = 8. The period of the moire pat¬ 
tern is about 100 times larger than the lattice constant 
of graphene. Here we fix the period of the potential, 
ao = 30 nm, and the amplitude of the potential, Vo = 20 
meV. We could then neglect the Landau level mixing, 
since it is very weak. The period is large enough to avoid 
the valley mixing. Hence, we could set the valley same 
as spin. The valley pseudo-spin is conserved in the HFA. 
In order to be consistent with the experimental results, 
we consider a £ [1,2], so 1/a £ [0.5,1] (a in Ref. 0 is 
1/a in this paper). Because of the fractal pattern, the 
noninteracting energy spectrum in the region a £ [1,2] 
is similar to that in the region a £ [0,1]. 

Finite-size study for v = — 1: For the finite-size study 
of the energy spectrum in the HFA, we fix the size of the 
sample, and change the magnetic field (or the parameter 
a). The size of the sample is 6ao x 6ao- For the filling 
factor v = — 1 in the Landau level TV = 0, there is AN^ 
degeneracy if we do not consider the Zeeman coupling, 
but there are only TV^ electrons in this LL, N s = N^. 
So in this sample, TV^ = [36/a], where [i] is the largest 
integer not exceeding x. 

In Fig. [0(a), we show the energy gaps with different a. 
The energy gap oscillates as the magnetic field increases. 
When a = 1.5 , which is equivalent to a = 0.5, a previous 
study [I| showed that the gap between the two bands 
would be open when the spin is polarized and one valley 
is half-filled in the Hartree approximation. In this work, 
we take the spin into consideration. The ground state is 
no longer spin polarized. The Zeeman coupling is very 



FIG. 1: The energy gaps for the filling factor v = — 1 with 
different 1/a. (a) The size of the sample is 6ao x 6ao. (b) An 
infinite sample. 



FIG. 2: Energy gaps (a) for v = 0 and (b) for v = 4. 


weak (about 1 meV) while the amplitude of the external 
potential is 20 meV. The potential is strong enough to 
mix different spins. The gap for a = 1.5 is very small in 
Fig.[H(a). This is because, intuitively, the two spins are 
mixed, and the corresponding four bands (two for each 
spin) in one valley are overlapped to close the Hartree 
gap. This mixed spin ground state will be discussed in 
detail in the infinite-size study below. 

Energy gaps in an infinite sample for v = — 1: The 
finite-size study however, may not be very reliable. In 
fact, solution of Eq. m is predicated on the size of 
the sample being infinite. In Fig. [T] (b), the energy gap 
also oscillates with 1/a. However, the amplitude and the 
peaks of the oscillation are changed a little. This might 
be because the system is size-dependent when the size is 
finite. For a = 1.5, the gap is also very small, which is 
simialr to that of the fintie-size calculation. Generally, 
the results of the two different calculations are similar. 
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FIG. 3: (Color online) The DOS for v = 4, and (a) for a = 1.1 
and (b) for a = 1.4. The arrows indicate the Fermi level 
locations. 


Experimental results (Fig. 4 (c) in 0 ) also show oscilla¬ 
tions in the energy gap, but the measured gap is nonzero 
for a = 2. 

Energy gaps in an infinite sample for v = 0,4: We 
calculate for filling factors v = 0,4, in the LLs N = 0,1, 
respectively. In these cases, each LL is half-filled, i.e., 
there are N s = 2 N$ electrons in the LL 0 or 1. The spin¬ 
less picture is obviously not satisfied in this case. We 
need to consider all the spins and valleys. Figure [2] (a) 
and Fig. [2] (b) show the energy gaps for the filling factors 
v = 0,4, respectively. These two curves are similar: the 
gaps are very small except at two points a = 1.4,1.6. 
Note that for a = 1.5, the energy gap is small, due to 
the same reason as what we explained for v = — 1. For 
the filling factor v = 0, the numerical results are different 
from the experimental results [l3|, where the energy gap 
curve looks like the energy cruve of the charged skyrmion 
excitation [24j . However, in our calculations we can not 
obtain such a skyrmion crystal ground state. It might be 
because the electron density is much higher than what 
the skyrmion crystal was found numerically [25j . The 
spin or pseudo-spin textures are suppressed by the high 
density electron gas: (pseudo-) spin flipping can not de¬ 
crease enough energy to create a (pseudo-) spin texture. 
The reason why our numerical results differ from the 
experiment is because perhaps the ground state in the 
N = 0 LL is not spin polarized without the external po¬ 
tential j3o}, and we only consider a spin polarized liquid 
ground state here. 

For v = 4, our numerical results are similar to those 
observed in the experiment [l3[. In the DOS, we clearly 
see the energy band structures. The DOS for a £ [1.5,2] 
is similar to the DOS for a = 3 — a, so that we neglect 
the DOS for a > 1.5. For simiplicity, we show the DOS 
curves for a = 1.1 and a = 1.4 in Figs. [3] (a) and [3] (b), 
respectively. For a = 1.1, there are ten bands, but the 
middle two bands touch at the Fermi level. The energy 



FIG. 4: (Color online) (a) Density profile for a = 1.4, v = 4, 
in the K valley. The density is in units of 1/(2tv£ 2 ). (b) The 
spin field, in units of ft/ (2n£ 2 ) , in the a direction Sk,z in the 
K valley for a = 1.4, v = 4 are showed. 


gap is almost zero. Note that some bands far away from 
the Fermi level split into two sub-bands, due to the Zee- 
man coupling. For a = 1.4, there is a gap between the 
middle two bands and all bands mix both spins, opening 
a gap (about 2.5 meV). 

We now define the spin field [in[ in valley p ( K or K') 

Sr],x + iSri,y = (P(?7, 4 , 0 , 4 ,) ( r )) > (14) 

Sr],z = (P(ri, T),(t),T) ( r ) ~ P(vM,(v,i) ( r )) ■ (15) 

The density in valley k is given by nk (r) = 

Y2 S (P{v,s),(ri,s) ( r )) ) where s is the spin index. The two 
valleys are completely equivalent in our numerical results, 
so only the order parameters in the K valley are shown 
in Fig. [4] The density profiles have the same geometry 
as the external potential. There is no valley coherence 
which is not showed in Fig. []J i.e. {p(k,s),(k' ,s')) = 0. 
The spin field contains no texture at all, S V}X = = 0, 

so the electron crystal is not a skyrmion crystal. Only 
the ^-components are nonzero, S VtZ ^ 0, and S v ^ z is also 
crystallized. 

Note that the maximum points of S v ^ z do not match 
the maximum points of the density, where the minimum 
points of the external potential are. At these points the 
potential decreases the kinetic energy for both spins. The 
electrons with both spins overcome the repulsive inter¬ 
action to be localized by the potential. The density of 
electrons is minimum at the sites where the energy of the 
potential is maximum. At the points where the density 
of electrons is maximum or minimum, the spin field S V}Z 
is minimum [the blue dots in Fig.[I](b)]. 

In conclusion, we have studied the interacting Hof- 
stadter’s butterfly states in a finite-size system in the 
HFA in order to study the spin/valley systems such as 
graphene. We also used a method where the sample is 
infinite. The energy gaps in the finite-size study agree 
with the results of the electron crystal qualitatively for 
filling factor v = — 1. The excitation gap oscillates with 
the increase of the magnetic field (or 1/a), similar to 
a recent experimental observation. For half fillings, we 
employ the crystal method to calculate the DOS of the 
system. In the n = 1 LL, the energy gap in a magnetic 
field (or 1/a) agrees well with the experimental results. 
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The osillation of the energy gap agrees qualitatively with 
that in the experiment, which the nonintearcing picture 
without spin or valley can not explain. Finally, we show 
the ground state of the electron gas in the Hofstadter’s 
butterfly state. The two spins are mixed while there is 
no valley coherence in the system. This crystal phase 
is neither like a Wigner crystal nor a skyrnrion crystal. 
The electron gas tends to be in a liquid phase, but the 


strong external potential crystallizes it. We propose that 
this method is able to study the interacting (in the HFA) 
Hofstadter’s butterfly states with different periodic po¬ 
tentials in an infinite system conveniently and efficiently, 
not only in graphene, but also in other Dirac materials. 

The work has been supported by the Canada Research 
Chairs Program of the Government of Canada. 
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